perm filename FITS.FAI[1,BGB] blob
sn#101481 filedate 1974-05-10 generic text, type C, neo UTF8
COMMENT ⊗ VALID 00007 PAGES
C REC PAGE DESCRIPTION
C00001 00001
C00002 00002 TITLE FITS - A LEAST SQUARES CUBIC FIT - NOVEMBER 1972.
C00004 00003 ACCUMULATE PRODUCTS OF THE COORDINATES OF THE DATA POINTS.
C00006 00004 Solve four equations in four unknowns by Gaussian Elimination.
C00008 00005 Find column-2 pivot for 3 by 3 matriX
C00009 00006 Find column-3 pivot for 2 by 2 matriX
C00010 00007 Back Substitution.
C00011 ENDMK
C⊗;
TITLE FITS - A LEAST SQUARES CUBIC FIT - NOVEMBER 1972.
COMMENT/
Y = (((B0*X + B1)*X + B2)*X + B3).
SY = B3*N + B2*SX + B1*X2 + B0*X3.
XY = B3*SX + B2*X2 + B1*X3 + B0*X4.
X2Y = B3*X2 + B2*X3 + B1*X4 + B0*X5.
X3Y = B3*X3 + B2*X4 + B1*X5 + B0*X6.
/
;COEFFICIENTS OF THE FOUR EQUATIONS.
ROW1: BLOCK 6
ROW2: BLOCK 6
ROW3: BLOCK 6
ROW4: BLOCK 6
;Use Accumulators for accumulating.
ACCUMULATORS{X,Y,SX,SY,X2,XY,X3,X2Y,X4,X3Y,X5,X6}
I←1
;CUBFIT(A,B,N)
; - halfword array of N points (SX,y) given in A.
; - returns four coefficients in real array B.
SUBR(CUBFIT)
dac 12,AC12
lac I,ARG3
lacn ARG1↔SOS
hrlm 0,I ;loop count and pointer.
lac[xwd 4,5]↔setz 4,↔blt 15 ;clear accumulators.
;ACCUMULATE PRODUCTS OF THE COORDINATES OF THE DATA POINTS.
L1: NIP X,(I) ↔ FSC X,233-12
NAP Y,(I) ↔ FSC Y,233-12
FADR SX,X
FADR SY,Y
FMPR Y,X ↔ FADR XY,Y
FMPR Y,X ↔ FADR X2Y,Y
FMPR Y,X ↔ FADR X3Y,Y
LAC X
FMPR X↔FADR X2,
FMPR X↔FADR X3,
FMPR X↔FADR X4,
FMPR X↔FADR X5,
FMPR X↔FADR X6,
L2: AOBJN I,L1
;SETUP THE COEFFICIENTS.
LAC 0,ARG1↔AOS↔FSC 0,233
MOVEI 1,ROW1
MOVEI 2,ROW2
DAC 0,(1)1↔DAC SX,(1)2↔DAC X2,(1)3↔DAC X3,(1)4↔DAC SY,(1)5
DAC SX,(2)1↔DAC X2,(2)2↔DAC X3,(2)3↔DAC X4,(2)4↔DAC XY,(2)5
MOVEI 3,ROW3
MOVEI 4,ROW4
DAC X2,(3)1↔DAC X3,(3)2↔DAC X4,(3)3↔DAC X5,(3)4↔DAC X2Y,(3)5
DAC X3,(4)1↔DAC X4,(4)2↔DAC X5,(4)3↔DAC X6,(4)4↔DAC X3Y,(4)5
;Solve four equations in four unknowns by Gaussian Elimination.
;TRIANGULARIZATION.
;Accumulators 1,2,3,4 are row pointers.
R←5 ;Reciprocal of the Pivot.
M←6 ;Multiplier of the row.
SETZM PARITY#
;Find column-1 pivot for 4 by 4 matrix.
T1: lacm (1)1↔lacm M,(2)1↔caml M↔jrst .+3↔exch 1,2↔setcmm PARITY
lacm (1)1↔lacm M,(3)1↔caml M↔jrst .+3↔exch 1,3↔setcmm PARITY
lacm (1)1↔lacm M,(4)1↔caml M↔jrst .+3↔exch 1,4↔setcmm PARITY
;Reciprocal of the Pivot.
movsi R,(1.0)↔fdvr R,(1)1
;Zero column-1 elements of rows 2, 3, & 4.
lac M,(2)1↔fmpr M,R↔ SETZM (2)1
lacn M↔fmpr (1)2↔fadrm (2)2
lacn M↔fmpr (1)3↔fadrm (2)3
lacn M↔fmpr (1)4↔fadrm (2)4
lacn M↔fmpr (1)5↔fadrm (2)5
lac M,(3)1↔fmpr M,R↔ SETZM (3)1
lacn M↔fmpr (1)2↔fadrm (3)2
lacn M↔fmpr (1)3↔fadrm (3)3
lacn M↔fmpr (1)4↔fadrm (3)4
lacn M↔fmpr (1)5↔fadrm (3)5
lac M,(4)1↔fmpr M,R↔ SETZM (4)1
lacn M↔fmpr (1)2↔fadrm (4)2
lacn M↔fmpr (1)3↔fadrm (4)3
lacn M↔fmpr (1)4↔fadrm (4)4
lacn M↔fmpr (1)5↔fadrm (4)5
;Normalize First Row.
FMPRM R,(1)1
fmprm R,(1)2
fmprm R,(1)3
fmprm R,(1)4
fmprm R,(1)5
;Find column-2 pivot for 3 by 3 matriX
T2: lacm(2)2↔lacm M,(3)2↔caml M↔jrst .+3↔exch 2,3↔setcmm PARITY
lacm(2)2↔lacm M,(4)2↔caml M↔jrst .+3↔exch 2,4↔setcmm PARITY
;Reciprocal of the pivot.
movsi R,(1.0)↔fdvr R,(2)2
;Zero column-2 elements of rows 3 & 4.
lac M,(3)2↔fmpr M,R ↔SETZM (3)2
lacn M↔fmpr (2)3↔fadrm (3)3
lacn M↔fmpr (2)4↔fadrm (3)4
lacn M↔fmpr (2)5↔fadrm (3)5
lac M,(4)2↔fmpr M,R ↔SETZM (4)2
lacn M↔fmpr (2)3↔fadrm (4)3
lacn M↔fmpr (2)4↔fadrm (4)4
lacn M↔fmpr (2)5↔fadrm (4)5
;Normalize Second Row.
FMPRM R,(2)2
fmprm R,(2)3
fmprm R,(2)4
fmprm R,(2)5
;Find column-3 pivot for 2 by 2 matriX
T3: lacm(3)3↔lacm M,(4)3↔caml M↔jrst .+3↔exch 3,4↔setcmm PARITY
;Reciprocal of the pivot.
movsi R,(1.0)↔fdvr R,(3)3
;Zero column-3 element of row 4.
lac M,(4)3↔fmpr M,R ↔SETZM (4)3
lacn M↔fmpr (3)4↔fadrm (4)4
lacn M↔fmpr (3)5↔fadrm (4)5
;Normalize Row-3.
FMPRM R,(3)3
fmprm R,(3)4
fmprm R,(3)5
;Find column-4 pivot for 1 by 1 matriX (trivail).
;Reciprocal of the pivot.
T4: movsi R,(1.0)↔fdvr R,(4)4
;Normalize Row-4.
FMPRM R,(4)4
fmprm R,(4)5
;Back Substitution.
T5: lacn (4)5↔fmpr (3)4↔fadrm (3)5
lacn (4)5↔fmpr (2)4↔fadrm (2)5
lacn (4)5↔fmpr (1)4↔fadrm (1)5
lacn (3)5↔fmpr (2)3↔fadrm (2)5
lacn (3)5↔fmpr (1)3↔fadrm (1)5
lacn (2)5↔fmpr (1)2↔fadrm (1)5
;Move results to Vector B.
lac 5,ARG2 ;B vector pointer.
skipe PARITY↔jrst T7
T6: lac(1)5↔dac(5)3
lac(2)5↔dac(5)2
lac(3)5↔dac(5)1
lac(4)5↔dac(5)0
lac 12,AC12↔pop3j
T7: lacn(1)5↔dac(5)3
lacn(2)5↔dac(5)2
lacn(3)5↔dac(5)1
lacn(4)5↔dac(5)0
lac 12,AC12↔pop3j
END